Chapter 5 Community composition
5.1 Taxonomy overview
5.1.1 Stacked barplot
# Merge data frames based on sample
merged_data<-genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
filter(count > 0) #filter 0 counts
# Create an interaction variable for time_point and sample
merged_data$interaction_var <- interaction(merged_data$sample, merged_data$time_point)
ggplot(merged_data, aes(x=interaction_var,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
scale_fill_manual(values=phylum_colors) +
facet_nested(. ~ type, scales="free") + #facet per day and treatment
guides(fill = guide_legend(ncol = 1)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.title.x = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
labs(fill="Phylum",y = "Relative abundance",x="Samples")+
scale_x_discrete(labels = function(x) gsub(".*\\.", "", x)) +
facet_wrap(~type, scales = "free", labeller = as_labeller(function(label) gsub(".*\\.", "", label))) #only show time_point label
5.1.2 Phylum relative abundances
phylum_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
left_join(genome_metadata, by = join_by(genome == genome)) %>%
group_by(sample,phylum) %>%
summarise(relabun=sum(count))
phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| phylum | mean | sd |
|---|---|---|
| p__Bacteroidota | 0.383073524 | 0.200556253 |
| p__Bacillota_A | 0.255324749 | 0.163209517 |
| p__Bacillota | 0.118782585 | 0.146832611 |
| p__Pseudomonadota | 0.094122159 | 0.156932597 |
| p__Campylobacterota | 0.053366123 | 0.093185234 |
| p__Verrucomicrobiota | 0.027449290 | 0.066220962 |
| p__Desulfobacterota | 0.023720834 | 0.036545754 |
| p__Chlamydiota | 0.010635637 | 0.059863092 |
| p__Fusobacteriota | 0.010544526 | 0.028384721 |
| p__Cyanobacteriota | 0.009261265 | 0.016525997 |
| p__Bacillota_C | 0.004741218 | 0.006655526 |
| p__Spirochaetota | 0.004033548 | 0.012340900 |
| p__Bacillota_B | 0.002480140 | 0.004938693 |
| p__Actinomycetota | 0.001243097 | 0.006365104 |
| p__Elusimicrobiota | 0.001221306 | 0.006520508 |
phylum_arrange <- phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun)) %>%
arrange(-mean) %>%
select(phylum) %>%
pull()
phylum_summary %>%
filter(phylum %in% phylum_arrange) %>%
mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
geom_jitter(alpha=0.5) +
theme_minimal() +
theme(legend.position="none") +
labs(y="Phylum",x="Relative abundance")
5.2 Taxonomy boxplot
5.2.1 Family
family_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,family) %>%
summarise(relabun=sum(count))
family_summary %>%
group_by(family) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| family | mean | sd |
|---|---|---|
| f__Bacteroidaceae | 2.228331e-01 | 0.1385585120 |
| f__Lachnospiraceae | 1.414562e-01 | 0.1083170968 |
| f__Tannerellaceae | 1.034197e-01 | 0.0797797418 |
| f__Helicobacteraceae | 5.290294e-02 | 0.0927590950 |
| f__Mycoplasmoidaceae | 3.716482e-02 | 0.0758142978 |
| f__Erysipelotrichaceae | 3.523239e-02 | 0.0451748821 |
| f__UBA3700 | 3.438944e-02 | 0.0558522623 |
| f__Marinifilaceae | 2.790560e-02 | 0.0271961340 |
| f__DTU072 | 2.716390e-02 | 0.0530196982 |
| f__Rikenellaceae | 2.712953e-02 | 0.0465396092 |
| f__Enterobacteriaceae | 2.707699e-02 | 0.0920895535 |
| f__Coprobacillaceae | 2.566962e-02 | 0.0894627842 |
| f__ | 2.527324e-02 | 0.0777603339 |
| f__Desulfovibrionaceae | 2.372083e-02 | 0.0365457537 |
| f__Ruminococcaceae | 1.871009e-02 | 0.0430525690 |
| f__UBA3830 | 1.588214e-02 | 0.0455209161 |
| f__Rhizobiaceae | 1.532665e-02 | 0.0768399323 |
| f__LL51 | 1.510265e-02 | 0.0608228261 |
| f__Akkermansiaceae | 1.234664e-02 | 0.0317643089 |
| f__Chlamydiaceae | 1.063564e-02 | 0.0598630919 |
| f__Fusobacteriaceae | 1.054453e-02 | 0.0283847213 |
| f__CAG-239 | 9.195743e-03 | 0.0151277476 |
| f__Enterococcaceae | 8.470645e-03 | 0.0467930117 |
| f__Gastranaerophilaceae | 7.774775e-03 | 0.0144711894 |
| f__Oscillospiraceae | 6.682761e-03 | 0.0078169630 |
| f__UBA1997 | 6.416374e-03 | 0.0310554945 |
| f__Streptococcaceae | 6.404337e-03 | 0.0343162283 |
| f__UBA1242 | 4.701177e-03 | 0.0154146977 |
| f__Brevinemataceae | 4.033548e-03 | 0.0123409004 |
| f__Acutalibacteraceae | 3.394636e-03 | 0.0109856431 |
| f__UBA660 | 3.215538e-03 | 0.0139953943 |
| f__Clostridiaceae | 2.859123e-03 | 0.0171837280 |
| f__RUG11792 | 2.834502e-03 | 0.0251869021 |
| f__Peptococcaceae | 2.480140e-03 | 0.0049386928 |
| f__MGBC116941 | 2.173027e-03 | 0.0094092981 |
| f__Acidaminococcaceae | 1.955375e-03 | 0.0050449819 |
| f__CAG-508 | 1.818740e-03 | 0.0064351937 |
| f__Anaerovoracaceae | 1.578874e-03 | 0.0036560185 |
| f__Moraxellaceae | 1.494257e-03 | 0.0097731073 |
| f__RUG14156 | 1.486490e-03 | 0.0045970591 |
| f__Staphylococcaceae | 1.369815e-03 | 0.0051004179 |
| f__Elusimicrobiaceae | 1.221306e-03 | 0.0065205081 |
| f__CAG-288 | 9.547359e-04 | 0.0060374248 |
| f__Anaerotignaceae | 9.043255e-04 | 0.0040584492 |
| f__CALVMC01 | 7.561589e-04 | 0.0043736226 |
| f__Eggerthellaceae | 6.446025e-04 | 0.0021324245 |
| f__Massilibacillaceae | 6.272186e-04 | 0.0016392219 |
| f__Mycobacteriaceae | 5.984944e-04 | 0.0060578862 |
| f__UBA1820 | 4.764220e-04 | 0.0013091229 |
| f__Arcobacteraceae | 4.631865e-04 | 0.0050473402 |
| f__CAG-274 | 4.546649e-04 | 0.0022091547 |
| f__Burkholderiaceae_C | 3.721451e-04 | 0.0048235514 |
| f__Muribaculaceae | 3.639429e-04 | 0.0009801614 |
| f__UBA932 | 3.439904e-04 | 0.0011614775 |
| f__Hepatoplasmataceae | 3.006899e-04 | 0.0038973865 |
| f__Rhodobacteraceae | 2.976706e-04 | 0.0038582514 |
| f__Weeksellaceae | 2.788124e-04 | 0.0031569776 |
| f__Eubacteriaceae | 1.656625e-04 | 0.0006747973 |
| f__Sphingobacteriaceae | 1.514738e-04 | 0.0012496721 |
| f__Devosiaceae | 1.498864e-04 | 0.0015139002 |
| f__Pumilibacteraceae | 1.285021e-04 | 0.0007668974 |
| f__WRAU01 | 9.660522e-05 | 0.0012521468 |
| f__Peptostreptococcaceae | 2.300953e-05 | 0.0002982376 |
family_arrange <- family_summary %>%
group_by(family) %>%
summarise(mean=sum(relabun)) %>%
arrange(-mean) %>%
select(family) %>%
pull()
family_summary %>%
left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
filter(family %in% family_arrange[1:20]) %>%
mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
scale_color_manual(values=phylum_colors[-8]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")
5.2.2 Genus
genus_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,genus) %>%
summarise(relabun=sum(count)) %>%
filter(genus != "g__")
genus_summary %>%
group_by(genus) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| genus | mean | sd |
|---|---|---|
| g__Bacteroides | 1.360888e-01 | 0.0920042645 |
| g__Parabacteroides | 9.738932e-02 | 0.0800669017 |
| g__Phocaeicola | 6.976407e-02 | 0.0795272001 |
| g__Mycoplasmoides | 3.069452e-02 | 0.0755136832 |
| g__Helicobacter_J | 3.026915e-02 | 0.0596129112 |
| g__Odoribacter | 2.567500e-02 | 0.0267638614 |
| g__Roseburia | 2.398913e-02 | 0.0560807126 |
| g__NHYM01 | 2.263379e-02 | 0.0799142802 |
| g__Alistipes | 2.221496e-02 | 0.0283922644 |
| g__Coprobacillus | 2.025846e-02 | 0.0881340464 |
| g__Agrobacterium | 1.532665e-02 | 0.0768399323 |
| g__Akkermansia | 1.234664e-02 | 0.0317643089 |
| g__Fusobacterium_A | 1.045087e-02 | 0.0283902792 |
| g__Kineothrix | 8.853496e-03 | 0.0410203841 |
| g__Proteus | 8.709520e-03 | 0.0683836665 |
| g__Dielma | 8.629168e-03 | 0.0090864942 |
| g__CAG-95 | 8.157620e-03 | 0.0205536177 |
| g__JAAYNV01 | 8.042189e-03 | 0.0195890923 |
| g__UBA866 | 7.303541e-03 | 0.0294323192 |
| g__Desulfovibrio | 7.059930e-03 | 0.0212042084 |
| g__Enterococcus | 7.042975e-03 | 0.0457933183 |
| g__Ureaplasma | 6.470299e-03 | 0.0138334733 |
| g__Lactococcus | 6.404337e-03 | 0.0343162283 |
| g__Lacrimispora | 6.100861e-03 | 0.0098371428 |
| g__Parabacteroides_B | 6.030362e-03 | 0.0100365790 |
| g__CALXRO01 | 5.800048e-03 | 0.0309414387 |
| g__Citrobacter | 5.751065e-03 | 0.0335518785 |
| g__NSJ-61 | 5.593207e-03 | 0.0199387071 |
| g__Breznakia | 5.537293e-03 | 0.0237686728 |
| g__Clostridium_AQ | 5.357893e-03 | 0.0121988801 |
| g__Salmonella | 5.164140e-03 | 0.0185065215 |
| g__Bilophila | 5.013506e-03 | 0.0088636099 |
| g__Hungatella_A | 4.840444e-03 | 0.0095762431 |
| g__MGBC136627 | 4.390777e-03 | 0.0163455561 |
| g__Escherichia | 4.213296e-03 | 0.0266876301 |
| g__UMGS1251 | 4.184603e-03 | 0.0072862642 |
| g__Hungatella | 4.140890e-03 | 0.0191331978 |
| g__Clostridium_Q | 4.035877e-03 | 0.0052190471 |
| g__Brevinema | 4.033548e-03 | 0.0123409004 |
| g__Thomasclavelia | 3.925809e-03 | 0.0109326053 |
| g__Scatousia | 3.671667e-03 | 0.0102995947 |
| g__Enterocloster | 3.641933e-03 | 0.0047397590 |
| g__Mailhella | 3.633580e-03 | 0.0102738868 |
| g__Copromonas | 3.620793e-03 | 0.0050252903 |
| g__Ventrimonas | 3.534268e-03 | 0.0071394041 |
| g__Caccovivens | 3.363217e-03 | 0.0122532776 |
| g__Lawsonia | 3.308592e-03 | 0.0117503555 |
| g__Fournierella | 3.236376e-03 | 0.0062446298 |
| g__Limenecus | 3.182108e-03 | 0.0065882600 |
| g__MGBC133411 | 3.060527e-03 | 0.0074762303 |
| g__Mucinivorans | 2.917358e-03 | 0.0374302862 |
| g__Sarcina | 2.859123e-03 | 0.0171837280 |
| g__Acetatifactor | 2.754437e-03 | 0.0058585690 |
| g__Eisenbergiella | 2.713468e-03 | 0.0068503944 |
| g__Bacteroides_G | 2.698691e-03 | 0.0347178966 |
| g__CAJLXD01 | 2.649496e-03 | 0.0087733544 |
| g__Blautia | 2.573938e-03 | 0.0061559204 |
| g__C-19 | 2.289060e-03 | 0.0048805592 |
| g__Butyricimonas | 2.230601e-03 | 0.0050201300 |
| g__Velocimicrobium | 2.214326e-03 | 0.0066941819 |
| g__CAZU01 | 2.209459e-03 | 0.0066119796 |
| g__MGBC116941 | 2.173027e-03 | 0.0094092981 |
| g__Intestinimonas | 2.094290e-03 | 0.0039467287 |
| g__Negativibacillus | 2.081393e-03 | 0.0055279017 |
| g__Rikenella | 1.997207e-03 | 0.0037146132 |
| g__Phascolarctobacterium | 1.955375e-03 | 0.0050449819 |
| g__RGIG6463 | 1.800497e-03 | 0.0039777189 |
| g__JALFVM01 | 1.752732e-03 | 0.0038715710 |
| g__Oscillibacter | 1.519449e-03 | 0.0025039951 |
| g__Pseudoflavonifractor | 1.511922e-03 | 0.0027764599 |
| g__Acinetobacter | 1.494257e-03 | 0.0097731073 |
| g__Citrobacter_A | 1.401214e-03 | 0.0060519635 |
| g__Staphylococcus | 1.369815e-03 | 0.0051004179 |
| g__RGIG4733 | 1.311551e-03 | 0.0040589196 |
| g__UBA1436 | 1.221306e-03 | 0.0065205081 |
| g__Lachnotalea | 1.215388e-03 | 0.0049307664 |
| g__Ruthenibacterium | 1.209176e-03 | 0.0053784424 |
| g__14-2 | 1.191694e-03 | 0.0096582871 |
| g__Beduini | 1.180938e-03 | 0.0025200906 |
| g__Scatocola | 1.127866e-03 | 0.0045101626 |
| g__Faecisoma | 1.092047e-03 | 0.0055756491 |
| g__Enterococcus_A | 1.090444e-03 | 0.0099443378 |
| g__Faecimonas | 9.939338e-04 | 0.0083324084 |
| g__RGIG9287 | 9.696200e-04 | 0.0093504118 |
| g__CAG-345 | 9.547359e-04 | 0.0060374248 |
| g__Blautia_A | 9.262837e-04 | 0.0029160490 |
| g__RGIG1896 | 8.393274e-04 | 0.0062956893 |
| g__CAG-269 | 8.029962e-04 | 0.0047313823 |
| g__Marseille-P3106 | 7.985520e-04 | 0.0017660762 |
| g__WRHT01 | 6.467835e-04 | 0.0027055875 |
| g__Eggerthella | 6.446025e-04 | 0.0021324245 |
| g__IOR16 | 6.437031e-04 | 0.0020707005 |
| g__Anaerotruncus | 6.302739e-04 | 0.0016992229 |
| g__RUG14156 | 6.188398e-04 | 0.0022208214 |
| g__CHH4-2 | 6.181620e-04 | 0.0020051727 |
| g__Corynebacterium | 5.984944e-04 | 0.0060578862 |
| g__Serratia_A | 5.895501e-04 | 0.0076414423 |
| g__CAG-56 | 4.944690e-04 | 0.0016386388 |
| g__Merdimorpha | 4.764220e-04 | 0.0013091229 |
| g__MGBC140009 | 4.707187e-04 | 0.0024136538 |
| g__CALURL01 | 4.662877e-04 | 0.0016783668 |
| g__Aliarcobacter | 4.631865e-04 | 0.0050473402 |
| g__Scatenecus | 4.547121e-04 | 0.0019816847 |
| g__RGIG8482 | 4.425249e-04 | 0.0029840978 |
| g__Enterobacter | 4.097683e-04 | 0.0041440045 |
| g__Klebsiella | 4.078573e-04 | 0.0049056069 |
| g__Caccenecus | 3.964657e-04 | 0.0017852972 |
| g__Alcaligenes | 3.721451e-04 | 0.0048235514 |
| g__Plesiomonas | 3.654875e-04 | 0.0027184625 |
| g__HGM05232 | 3.639429e-04 | 0.0009801614 |
| g__JAHHSE01 | 3.613502e-04 | 0.0014942841 |
| g__Egerieousia | 3.439904e-04 | 0.0011614775 |
| g__Emergencia | 3.433949e-04 | 0.0017500504 |
| g__Enterococcus_B | 3.372270e-04 | 0.0022331962 |
| g__Stoquefichus | 3.044085e-04 | 0.0020563925 |
| g__Hepatoplasma | 3.006899e-04 | 0.0038973865 |
| g__Paracoccus | 2.976706e-04 | 0.0038582514 |
| g__Moheibacter | 2.788124e-04 | 0.0031569776 |
| g__Scatomorpha | 2.656736e-04 | 0.0010212728 |
| g__UBA7185 | 2.448818e-04 | 0.0014600491 |
| g__Eubacterium | 1.656625e-04 | 0.0006747973 |
| g__Sphingobacterium | 1.514738e-04 | 0.0012496721 |
| g__Devosia | 1.498864e-04 | 0.0015139002 |
| g__Anaerosporobacter | 1.462768e-04 | 0.0012784872 |
| g__Caccomorpha | 1.391355e-04 | 0.0010571570 |
| g__UBA2658 | 1.315233e-04 | 0.0007225792 |
| g__Protoclostridium | 1.285021e-04 | 0.0007668974 |
| g__Angelakisella | 1.276029e-04 | 0.0009248545 |
| g__Cetobacterium_A | 9.365634e-05 | 0.0008744341 |
| g__Rahnella | 6.509221e-05 | 0.0008436915 |
| g__Peptostreptococcus | 2.300953e-05 | 0.0002982376 |
genus_arrange <- genus_summary %>%
group_by(genus) %>%
summarise(mean=sum(relabun)) %>%
filter(genus != "g__")%>%
arrange(-mean) %>%
select(genus) %>%
mutate(genus= sub("^g__", "", genus)) %>%
pull()
genus_summary %>%
left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
mutate(genus= sub("^g__", "", genus)) %>%
filter(genus %in% genus_arrange[1:20]) %>%
mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")
5.3 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
functional <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, dist = dist) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(functional = 1) %>%
rownames_to_column(var = "sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample)) %>%
full_join(functional, by = join_by(sample == sample))5.3.1 Wild samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric ) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
5.3.2 Acclimation samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
5.3.3 Antibiotics samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="2_Antibiotics") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
5.3.4 Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="3_Transplant1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
5.3.5 Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="4_Transplant2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
5.3.6 Post-Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
5.3.7 Post-Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
5.4 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, tree = genome_tree)
beta_q1f <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, dist = dist)5.5 Permanovas
5.5.1 1. Are the wild populations similar?
5.5.1.1 Wild: P.muralis vs P.liolepis
wild <- meta %>%
filter(time_point == "0_Wild")
# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))
wild_nmds <- sample_metadata %>%
filter(time_point == "0_Wild")5.5.1.3 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.002848 0.0028483 0.2723 999 0.624
Residuals 26 0.272017 0.0104622
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.619
Podarcis_muralis 0.60624 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.631953 | 0.207198 | 6.795072 | 0.001 |
| Residual | 26 | 6.244347 | 0.792802 | NA | NA |
| Total | 27 | 7.876300 | 1.000000 | NA | NA |
5.5.1.4 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.002893 0.0028935 0.2628 999 0.598
Residuals 26 0.286298 0.0110115
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.593
Podarcis_muralis 0.61255 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 2.018797 | 0.2566342 | 8.976049 | 0.001 |
| Residual | 26 | 5.847641 | 0.7433658 | NA | NA |
| Total | 27 | 7.866438 | 1.0000000 | NA | NA |
5.5.1.5 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.07936 0.079362 4.7909 999 0.037 *
Residuals 26 0.43070 0.016565
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.033
Podarcis_muralis 0.037786 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.3786052 | 0.2108638 | 6.947419 | 0.001 |
| Residual | 26 | 1.4168908 | 0.7891362 | NA | NA |
| Total | 27 | 1.7954960 | 1.0000000 | NA | NA |
5.5.1.6 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.016513 0.016513 1.4504 999 0.254
Residuals 26 0.296006 0.011385
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.277
Podarcis_muralis 0.23931 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0787916 | 0.1594096 | 4.930642 | 0.072 |
| Residual | 26 | 0.4154800 | 0.8405904 | NA | NA |
| Total | 27 | 0.4942716 | 1.0000000 | NA | NA |
beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
5.5.2 2. Do the antibiotics work?
5.5.2.1 Acclimation vs antibiotics
treat <- meta %>% #antibiotics samples giving error when plotting
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))
treat_nmds <- sample_metadata %>%
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")5.5.2.2 Number of samples used
[1] 52beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)5.5.2.2.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.027648 0.0276478 6.8746 999 0.01 **
Residuals 50 0.201088 0.0040218
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.015
2_Antibiotics 0.011555 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 1.952865 | 0.0941315 | 6.574358 | 0.001 |
| species | 1 | 2.299135 | 0.1108223 | 7.740081 | 0.001 |
| individual | 26 | 9.662166 | 0.4657330 | 1.251072 | 0.003 |
| Residual | 23 | 6.831983 | 0.3293133 | NA | NA |
| Total | 51 | 20.746150 | 1.0000000 | NA | NA |
5.5.2.2.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.043782 0.043782 8.3826 999 0.004 **
Residuals 50 0.261147 0.005223
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.005
2_Antibiotics 0.0056034 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.078576 | 0.1048012 | 9.380227 | 0.001 |
| species | 1 | 3.046369 | 0.1535971 | 13.747699 | 0.001 |
| individual | 26 | 9.611970 | 0.4846327 | 1.668348 | 0.001 |
| Residual | 23 | 5.096598 | 0.2569690 | NA | NA |
| Total | 51 | 19.833512 | 1.0000000 | NA | NA |
5.5.2.2.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.60296 0.60296 37.862 999 0.001 ***
Residuals 50 0.79626 0.01593
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.001
2_Antibiotics 1.2645e-07 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 1.8994052 | 0.2145480 | 20.282274 | 0.001 |
| species | 1 | 0.8219541 | 0.0928441 | 8.777010 | 0.001 |
| individual | 26 | 3.9777800 | 0.4493115 | 1.633678 | 0.003 |
| Residual | 23 | 2.1539163 | 0.2432964 | NA | NA |
| Total | 51 | 8.8530557 | 1.0000000 | NA | NA |
5.5.2.2.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.18583 0.185833 5.3501 999 0.03 *
Residuals 50 1.73673 0.034735
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.032
2_Antibiotics 0.024874 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 1.935980 | 0.3884761 | 38.1894171 | 0.001 |
| species | 1 | 0.003156 | 0.0006333 | 0.0622557 | 0.809 |
| individual | 26 | 1.878423 | 0.3769266 | 1.4251551 | 0.172 |
| Residual | 23 | 1.165966 | 0.2339640 | NA | NA |
| Total | 51 | 4.983525 | 1.0000000 | NA | NA |
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_treat <- beta_div_func_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
5.5.3 3. Do the FMT work?
5.5.3.1 Comparison between FMT2 vs Post-FMT2
transplant3 <- meta %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")
transplant3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts)),sort(as.character(rownames(transplant3))))
transplant3_nmds <- sample_metadata %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")5.5.3.2 Number of samples used
[1] 53beta_div_richness_transplant3<-hillpair(data=transplant3.counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3.counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3.counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3.counts, q=1, dist=dist)5.5.3.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.9809525 | 0.0778442 | 5.775712 | 0.001 |
| time_point | 1 | 0.7092107 | 0.0562799 | 4.175734 | 0.001 |
| type | 1 | 1.4479860 | 0.1149060 | 8.525540 | 0.001 |
| individual | 25 | 6.9157236 | 0.5488022 | 1.628753 | 0.001 |
| Residual | 15 | 2.5476146 | 0.2021678 | NA | NA |
| Total | 43 | 12.6014875 | 1.0000000 | NA | NA |
5.5.3.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.1101711 | 0.0895457 | 7.879642 | 0.001 |
| time_point | 1 | 0.7363935 | 0.0593971 | 5.226687 | 0.001 |
| type | 1 | 1.9027421 | 0.1534740 | 13.505059 | 0.001 |
| individual | 25 | 6.5351386 | 0.5271203 | 1.855374 | 0.001 |
| Residual | 15 | 2.1133659 | 0.1704628 | NA | NA |
| Total | 43 | 12.3978112 | 1.0000000 | NA | NA |
5.5.3.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.1341431 | 0.1005265 | 6.921553 | 0.001 |
| time_point | 1 | 0.1159136 | 0.0868654 | 5.980943 | 0.001 |
| type | 1 | 0.1423573 | 0.1066822 | 7.345390 | 0.001 |
| individual | 25 | 0.6512838 | 0.4880704 | 1.344205 | 0.067 |
| Residual | 15 | 0.2907075 | 0.2178554 | NA | NA |
| Total | 43 | 1.3344053 | 1.0000000 | NA | NA |
5.5.3.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | -0.0031096 | -0.0030953 | -0.1177247 | 0.824 |
| time_point | 1 | -0.0045864 | -0.0045653 | -0.1736359 | 0.833 |
| type | 1 | 0.0775554 | 0.0771987 | 2.9361453 | 0.141 |
| individual | 25 | 0.5385511 | 0.5360739 | 0.8155530 | 0.621 |
| Residual | 15 | 0.3962105 | 0.3943880 | NA | NA |
| Total | 43 | 1.0046211 | 1.0000000 | NA | NA |
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
p1<-beta_neutral_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylo_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_func_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")
5.5.4 4. Are there differences between the control and the treatment group?
5.5.4.2 Number of samples used
[1] 27beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)5.5.4.2.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.015090 0.0075452 2.1059 999 0.136
Residuals 24 0.085988 0.0035828
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.004000 0.638
Hot_control 0.009502 0.220
Treatment 0.624826 0.248542 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.6488906 | 0.0757032 | 2.115638 | 0.002 |
| type | 1 | 0.5615418 | 0.0655126 | 1.830847 | 0.008 |
| Residual | 24 | 7.3610760 | 0.8587842 | NA | NA |
| Total | 26 | 8.5715084 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.3294975 0.9684059 0.06064512 0.499 1.000
2 Control vs Hot_control 1 0.7123847 2.3071902 0.11949901 0.003 0.009 *
3 Treatment vs Hot_control 1 0.4044284 1.3387544 0.07721168 0.054 0.162 5.5.4.2.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.006755 0.0033773 0.4095 999 0.69
Residuals 24 0.197956 0.0082482
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.30200 0.957
Hot_control 0.30023 0.565
Treatment 0.95485 0.53193 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.7984743 | 0.104299 | 3.065171 | 0.001 |
| type | 1 | 0.6051778 | 0.079050 | 2.323148 | 0.006 |
| Residual | 24 | 6.2519772 | 0.816651 | NA | NA |
| Total | 26 | 7.6556293 | 1.000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.2860314 0.9858296 0.06166897 0.427 1.000
2 Control vs Hot_control 1 0.8705308 3.2919398 0.16222894 0.003 0.009 *
3 Treatment vs Hot_control 1 0.4377241 1.6307998 0.09249721 0.037 0.111 5.5.4.2.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.00404 0.0020184 0.1345 999 0.909
Residuals 24 0.36024 0.0150100
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.84200 0.693
Hot_control 0.81332 0.812
Treatment 0.63160 0.76897 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0700374 | 0.0635900 | 1.6594454 | 0.15 |
| type | 1 | 0.0184254 | 0.0167292 | 0.4365646 | 0.80 |
| Residual | 24 | 1.0129278 | 0.9196808 | NA | NA |
| Total | 26 | 1.1013906 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.01869516 0.4206549 0.02727867 0.781 1.000
2 Control vs Hot_control 1 0.07762661 2.2806818 0.11828844 0.036 0.108
3 Treatment vs Hot_control 1 0.03350315 0.6872008 0.04118131 0.668 1.000 5.5.4.2.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.01287 0.0064374 0.4319 999 0.66
Residuals 24 0.35771 0.0149044
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.38600 0.758
Hot_control 0.38202 0.598
Treatment 0.75152 0.59510 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | -0.0048931 | -0.0058792 | -0.1628739 | 0.826 |
| type | 1 | 0.1161466 | 0.1395538 | 3.8660898 | 0.100 |
| Residual | 24 | 0.7210172 | 0.8663254 | NA | NA |
| Total | 26 | 0.8322707 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.14708071 6.531088 0.30333294 0.043 0.129
2 Control vs Hot_control 1 0.03067531 1.027084 0.05697451 0.341 1.000
3 Treatment vs Hot_control 1 0.04199233 1.256701 0.07282392 0.301 0.903 beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
p1<-beta_neutral_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")
5.5.4.4 Number of samples used
[1] 28beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)5.5.4.4.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.001996 0.0009978 0.2043 999 0.825
Residuals 25 0.122093 0.0048837
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.70400 0.805
Hot_control 0.67713 0.609
Treatment 0.79246 0.58926 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | 1.542992 | 0.1948479 | 3.025016 | 0.001 |
| Residual | 25 | 6.375966 | 0.8051521 | NA | NA |
| Total | 27 | 7.918958 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.2206497 0.7906964 0.04709134 0.862 1.000
2 Control vs Hot_control 1 0.7211095 2.7164967 0.13777786 0.001 0.003 *
3 Treatment vs Hot_control 1 0.7163430 2.6326288 0.13409456 0.001 0.003 *5.5.4.4.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.007233 0.0036166 0.6989 999 0.497
Residuals 25 0.129365 0.0051746
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.51500 0.645
Hot_control 0.49241 0.283
Treatment 0.65999 0.28046 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | 1.962645 | 0.2556667 | 4.293552 | 0.001 |
| Residual | 25 | 5.713932 | 0.7443333 | NA | NA |
| Total | 27 | 7.676577 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.1081320 0.3901033 0.02380115 0.998 1.000
2 Control vs Hot_control 1 0.7336222 2.8624966 0.14411565 0.001 0.003 *
3 Treatment vs Hot_control 1 0.7084005 2.6970385 0.13692609 0.001 0.003 *5.5.4.4.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.000406 0.0002032 0.0494 999 0.963
Residuals 25 0.102865 0.0041146
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.94000 0.831
Hot_control 0.93747 0.787
Treatment 0.84169 0.75551 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | 0.1793566 | 0.2173098 | 3.470559 | 0.001 |
| Residual | 25 | 0.6459931 | 0.7826902 | NA | NA |
| Total | 27 | 0.8253496 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.005034708 0.178068 0.01100675 0.979 1.000
2 Control vs Hot_control 1 0.099703372 3.519604 0.17152397 0.003 0.009 *
3 Treatment vs Hot_control 1 0.079909458 2.912000 0.14624347 0.004 0.012 .5.5.4.4.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.00836 0.0041823 0.2225 999 0.844
Residuals 25 0.46991 0.0187964
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.58600 0.679
Hot_control 0.53209 0.855
Treatment 0.60231 0.82567 | Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | 0.0031465 | 0.0045388 | 0.0569937 | 0.88 |
| Residual | 25 | 0.6900904 | 0.9954612 | NA | NA |
| Total | 27 | 0.6932368 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.010979753 0.6382789 0.038362075 0.424 1
2 Control vs Hot_control 1 0.003756213 0.1351199 0.007885556 0.702 1
3 Treatment vs Hot_control 1 0.019226092 0.5508465 0.031385752 0.497 1 beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="top")
p1<-beta_neutral_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")